version 18.0               // version control
set processors 8           // to ensure replicability across different numbers of cores
clear all                  // clear existing data
macro drop _all            // and macros, clean slate
set seed 20220909          // set seed
/* 
Notes: 
1) Must put "set sortseed #" right on top of EACH "telasso" (rather than at the top of the do file) to ensure the same number of selected controls will be reported.
2) For detailed explanations, see Stata's technical support's email chain on 10/10/2023. 
*/

local pgm  "dst-Appendix_Table_A5a_depression_telasso_alt_spec"      // file name
local who  "Muzhe Yang"                                              // author
local dte  "2025-01-22"                                              // created date
local dte2 "`c(current_date)'"                                       // last run date
local tag  "`pgm'.do, created by `who' on `dte', last run on `dte2'"

capture log close
log using "code\analysis\tables\appendix\\`pgm'.txt", text replace 
display "`tag'"

**# data prep ------------------------------------------------------------------------------------

use "data_clean\dst-data04_for_estimation_within_250_miles", clear
summ dist_to_border
tab wave, missing 
des depression
keep if wave == 2
codebook TractFIPS

*## (1) create the 9 cells, following page 215 of the following paper:
*## Giuntella, O. and F. Mazzonna (2019). "Sunset Time and the Economic Effects of Social Jetlag: Evidence from US Time Zone Borders." Journal of Health Economics 65: 210-226.

assert !missing(centroid_lat)
assert !missing(region)
gen     cell = 1 if centroid_lat > 40                        & region == "eastern and central"
replace cell = 2 if (34 < centroid_lat & centroid_lat <= 40) & region == "eastern and central"
replace cell = 3 if centroid_lat <= 34                       & region == "eastern and central"
replace cell = 4 if centroid_lat > 40                        & region == "central and mountain"
replace cell = 5 if (34 < centroid_lat & centroid_lat <= 40) & region == "central and mountain"
replace cell = 6 if centroid_lat <= 34                       & region == "central and mountain"
replace cell = 7 if centroid_lat > 40                        & region == "mountain and pacific"
replace cell = 8 if (34 < centroid_lat & centroid_lat <= 40) & region == "mountain and pacific"
replace cell = 9 if centroid_lat <= 34                       & region == "mountain and pacific"
tab cell time_zone, missing
summ centroid_lat if cell == 1 | cell == 4 | cell == 7
summ centroid_lat if cell == 2 | cell == 5 | cell == 8
summ centroid_lat if cell == 3 | cell == 6 | cell == 9

*## (2) control variables

global x "pop white_pct black_pct hispanic_pct pop_under_18yr_pct pop_65yr_over_pct age_median educ_hs_pct educ_coll_pct married_pct hh_size hh_income_median home_value_median no_health_ins_pct unemploy_pct"
global cvars    c.($x dist_to_border centroid_lat daylight_dur_max)##c.($x dist_to_border centroid_lat daylight_dur_max)
global controls $cvars i.cell i.cell#c.($x dist_to_border centroid_lat daylight_dur_max) 

des  $x dist_to_border centroid_lat daylight_dur_max
summ $x dist_to_border centroid_lat daylight_dur_max

**# estimation prep -------------------------------------------------

encode CountyFIPS, gen(county_fips)
quietly tab cell, gen(cell_)
egen nomiss = rowmiss($x dist_to_border centroid_lat daylight_dur_max)
assert !missing(StateAbbr)
local radius = 50
gen sample_selection = (nomiss == 0 & StateAbbr != "AZ" & dist_to_border <= `radius')

**# table ---------------------------------------------------------------------------------

set sortseed 12102023
telasso (depression $controls) (treat $controls) if sample_selection == 1, vce(cluster county_fips) 
scalar k_controls = e(k_controls)
scalar k_controls_sel = e(k_controls_sel)
scalar method = "Plugin"
lassoinfo 
scalar lambda_y0 = r(table)["treat = 0", "lambda"]
scalar lambda_y1 = r(table)["treat = 1", "lambda"]
scalar lambda_treat = r(table)["treat", "lambda"]
scalar lati_sel = "LASSO"
scalar demo_sel = "LASSO"
scalar dist_sel = "LASSO"
scalar dayl_sel = "LASSO"
scalar cell_sel = "LASSO"
quietly etable, replace ///
        column(index) ///
		keep(r1vs0.treat) ///
        cstat(_r_b,  nformat(%9.3f)) ///
		cstat(_r_se, nformat(%9.3f)) ///
		mstat(N, nformat(%9.0fc) label("Number of observations")) ///
		mstat(k_controls, nformat(%9.0fc) label("Number of potential predictor variables")) ///
		mstat(k_controls_sel, nformat(%9.0fc) label("Number of selected predictor variables")) ///
		mstat(case_method = method, label("Method for selecting the value of the LASSO penalty parameter (i.e., lambda)")) ///
		mstat(case_lambda_y0 = lambda_y0, nformat(%9.3f) label("Lambda for the untreated outcome model")) ///
		mstat(case_lambda_y1 = lambda_y1, nformat(%9.3f) label("Lambda for the treated outcome model")) ///
		mstat(case_lambda_treat = lambda_treat, nformat(%9.3f) label("Lambda for the treatment model")) ///
		mstat(case_lati_sel = lati_sel, label("Latitude of the centroid of a census tract, selection done by LASSO or always included")) ///
		mstat(case_demo_sel = demo_sel, label("Demographic variables, selection done by LASSO or always included")) ///
		mstat(case_dist_sel = dist_sel, label("Distance to the time zone border, selection done by LASSO or always included")) ///
		mstat(case_dayl_sel = dayl_sel, label("Maximum of daylight hours per day within a year at the county level, selection done by LASSO or always included")) ///
		mstat(case_cell_sel = cell_sel, label("Cell dummy variables, selection done by LASSO or always included")) ///
        stars(0.10 "*" 0.05 "**" 0.01 "***", attach(_r_b) decreasing pvname("p-value") nformat(%9.2f)) showstars showstarsnote ///
        novarlabel nocenter 

set sortseed 12102023
telasso (depression $controls, lasso(selection(bic))) (treat $controls, lasso(selection(bic))) if sample_selection == 1, vce(cluster county_fips) 
scalar k_controls = e(k_controls)
scalar k_controls_sel = e(k_controls_sel)
scalar method = "BIC"
lassoinfo 
scalar lambda_y0 = r(table)["treat = 0", "lambda"]
scalar lambda_y1 = r(table)["treat = 1", "lambda"]
scalar lambda_treat = r(table)["treat", "lambda"]
scalar lati_sel = "LASSO"
scalar demo_sel = "LASSO"
scalar dist_sel = "LASSO"
scalar dayl_sel = "LASSO"
scalar cell_sel = "LASSO"
quietly etable, append

set sortseed 12102023
telasso (depression $controls, lasso(selection(cv))) (treat $controls, lasso(selection(bic))) if sample_selection == 1, vce(cluster county_fips) rseed(10101) 
scalar k_controls = e(k_controls)
scalar k_controls_sel = e(k_controls_sel)
scalar method = "CV"
lassoinfo 
scalar lambda_y0 = r(table)["treat = 0", "lambda"]
scalar lambda_y1 = r(table)["treat = 1", "lambda"]
scalar lambda_treat = r(table)["treat", "lambda"]
scalar lati_sel = "LASSO"
scalar demo_sel = "LASSO"
scalar dist_sel = "LASSO"
scalar dayl_sel = "LASSO"
scalar cell_sel = "LASSO"
quietly etable, append

set sortseed 12102023
telasso (depression $controls, lasso(selection(adaptive))) (treat $controls, lasso(selection(bic))) if sample_selection == 1, vce(cluster county_fips) rseed(10101) 
scalar k_controls = e(k_controls)
scalar k_controls_sel = e(k_controls_sel)
scalar method = "Adaptive CV"
lassoinfo 
scalar lambda_y0 = r(table)["treat = 0", "lambda"]
scalar lambda_y1 = r(table)["treat = 1", "lambda"]
scalar lambda_treat = r(table)["treat", "lambda"]
scalar lati_sel = "LASSO"
scalar demo_sel = "LASSO"
scalar dist_sel = "LASSO"
scalar dayl_sel = "LASSO"
scalar cell_sel = "LASSO"
quietly etable, append

set sortseed 12102023
telasso (depression $controls, ainclude($x dist_to_border centroid_lat daylight_dur_max cell_1-cell_9)) (treat $controls, ainclude(dist_to_border centroid_lat daylight_dur_max)) if sample_selection == 1, vce(cluster county_fips)
scalar k_controls = e(k_controls)
scalar k_controls_sel = e(k_controls_sel)
scalar method = "Plugin"
lassoinfo 
scalar lambda_y0 = r(table)["treat = 0", "lambda"]
scalar lambda_y1 = r(table)["treat = 1", "lambda"]
scalar lambda_treat = r(table)["treat", "lambda"]
scalar lati_sel = "Always included in the outcome and the treatment models"
scalar demo_sel = "Always included in the outcome model"
scalar dist_sel = "Always included in the outcome and the treatment models"
scalar dayl_sel = "Always included in the outcome and the treatment models"
scalar cell_sel = "Always included in the outcome model"
quietly etable, append

collect title "Table: The Impact of Circadian Misalignment on Depression Estimated by the ML with Alternative Specifications" 

collect layout 
collect addtags top_row[1]
collect layout (colname[r1vs0.treat]#result[_r_b _r_se] result[N k_controls k_controls_sel ///
                                                               case_method case_lambda_y0 case_lambda_y1 case_lambda_treat case_lati_sel case_demo_sel case_dist_sel case_dayl_sel case_cell_sel]) ///
               (top_row[1]#cmdset#stars)

collect label levels top_row ///
        1 "Treat = 1 for census tracts located east of the time zone border; Treat = 0 for census tracts located west of the time zone border", modify
collect label levels cmdset 1 "(1)" 2 "(2)" 3 "(3)" 4 "(4)" 5 "(5)", modify
collect label levels colname r1vs0.treat "Treat (1/0)", modify 

collect style header top_row, level(label)
collect style header cmdset, level(label)
collect style header colname, level(label)
collect style cell result[_r_b dep_mean N k_controls k_controls_sel case_lambda_y0 case_lambda_y1 case_lambda_treat], sformat(%s)
collect style cell result[effect_relative_magnitude], sformat("%s%%")
collect style cell border_block[corner], border(top, pattern(double))
collect style cell border_block[column-header], border(top, pattern(double))
collect style cell border_block[row-header], border(bottom, pattern(double))
collect style cell border_block[item], border(bottom, pattern(double))

collect export "code\analysis\tables\appendix\\`pgm'.xlsx", replace

log close
exit